1 Introduction

First, we again create 1 complete dataset, which is done by means of imputing the boys dataset in mice (Van Buuren & Groothuis-Oudshoorn, 2011) with default settings, and complete the dataset. Then, based on this one true dataset, we extract the regression coefficients of a linear regression model in which weight is regressed on age and height, that are used to compare the estimates of the synthetic versions of the datasets. The true regression coefficients are equal to \(\beta_0 = -9.216\), \(\beta_{age} = 2.327\), \(\beta_{height} = 0.191\). Then, for a total of 200 iterations, we impute a matrix with the same dimensions as the boys dataset (that is, the sample size \(n = 748\) and the number of predictors \(k = 9\)). This is done once with the default settings, that is, with pmm for the continuous variables, polr for binary variables and polyreg for ordinal variables with more than two categories. Furthermore, bmi is imputed passively, since it consists of a fixed combination of height and weight bmi = (wgt / (hgt/100)^2. Furthermore, the predictor matrix is adjusted so that imputations for bmi do not flow back into predictions of hgt and wgt. Then, for every synthetic dataset, the estimates, and corresponding variance of the estimates are calculated, as well as the \(95\%~CI\) are calculated. Furthermore, an indicator is added, whether or not the confidence interval includes the true parameter estimates. Note that these confidence intervals are still based on imputations for missing data, instead of imputations for synthetic data. This is due to the fact that the calculation for the degrees of freedom for the estimates based imputing synthetic values might yield degrees of freedom that are so small that the corresponding \(95\% ~ CI\) yields boundary values of \(-\infty\) and \(\infty\).


2 Simulations


2.1 Single dataset approach


2.1.1 True results

broom::tidy(truemodel, conf.int=TRUE) %>%
  mutate(CIW = conf.high - conf.low) %>% 
  kable(digits = 3,
        caption = "Results of a single linear regression model on the completed dataset.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Results of a single linear regression model on the completed dataset.
term estimate std.error statistic p.value conf.low conf.high CIW
(Intercept) -9.216 2.033 -4.534 0 -13.207 -5.226 7.982
age 2.327 0.189 12.335 0 1.956 2.697 0.741
hgt 0.191 0.028 6.843 0 0.136 0.246 0.110

2.1.2 Use mice to synthesize and to pool


2.1.2.1 Pool the results of the synthetic datasets with the default mice settings.

default_fit <- syns_def %>%
  map(function(x) x %$% lm(wgt ~ age + hgt) %>% pool %>% summary) %>%
  map_dfr(function(x) {
    var      <- x[,1]
    true_est <- coefs
    est      <- x[,2]
    true_se  <- sqrt(diag(vcov(truemodel)))
    se       <- x[,3]
    df       <- x[,5]
    lower    <- est + se * qt(.025, df)
    upper    <- est + se * qt(.975, df)
    cov      <- lower < coefs & coefs < upper
  
    bind_cols(var = var, true_est = true_est, est = est, true_se = true_se, 
                               se = se, df = df, lower = lower, upper = upper, cov = cov)
    })

results_def <- default_fit %>%
  group_by(var) %>%
  summarise("True Est" = unique(true_est),
            "Syn Est"  = mean(est),
            "Bias"     = mean(est - true_est),
            "True SE"  = unique(true_se),
            "Syn SE"   = mean(se),
            "df"       = mean(df),
            "Lower"    = mean(lower),
            "Upper"    = mean(upper),
            "CIW"      = mean(upper - lower),
            "Coverage" = mean(cov))

2.1.2.2 Pool the results of the synthetic dataset with CART.

cart_fit <- syns_cart %>%
  map(function(x) x %$% lm(wgt ~ age + hgt) %>% pool %>% summary) %>%
  map_dfr(function(x) {
    var      <- x[,1]
    true_est <- coefs
    est      <- x[,2]
    true_se  <- sqrt(diag(vcov(truemodel)))
    se       <- x[,3]
    df       <- x[,5]
    lower    <- est + se * qt(.025, df)
    upper    <- est + se * qt(.975, df)
    cov      <- lower < coefs & coefs < upper

    bind_cols(var = var, true_est = true_est, est = est, true_se = true_se,
              se = se, df = df, lower = lower, upper = upper, cov = cov)
    })

 results_cart <- cart_fit %>%
   group_by(var) %>%
   summarise("True Est" = unique(true_est),
             "Syn Est"  = mean(est),
             "Bias"     = mean(est - true_est),
             "True SE"  = unique(true_se),
             "Syn SE"   = mean(se),
             "df"       = mean(df),
             "Lower"    = mean(lower),
             "Upper"    = mean(upper),
             "CIW"      = mean(upper - lower),
             "Coverage" = mean(cov))

2.1.2.3 Results for both synthesization methods

bind_rows("Mice default" = results_def,
          "Mice with cart" = results_cart, .id = "Imputation method") %>%
  kable(digits = 3,
        caption = "Results of with mice default settings, and with CART settings, for the imputation of synthetic data, with the in-build mice pooling function.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Results of with mice default settings, and with CART settings, for the imputation of synthetic data, with the in-build mice pooling function.
Imputation method var True Est Syn Est Bias True SE Syn SE df Lower Upper CIW Coverage
Mice default (Intercept) -9.216 -2.160 7.057 2.033 2.941 28.084 -8.657 4.338 12.995 0.338
Mice default age 2.327 2.948 0.621 0.189 0.285 18.624 2.304 3.591 1.287 0.490
Mice default hgt 0.191 0.094 -0.097 0.028 0.041 22.553 0.003 0.186 0.183 0.388
Mice with cart (Intercept) -9.216 -7.602 1.614 2.033 2.664 50.879 -13.212 -1.992 11.220 1.000
Mice with cart age 2.327 2.463 0.136 0.189 0.256 41.663 1.917 3.008 1.091 1.000
Mice with cart hgt 0.191 0.169 -0.022 0.028 0.037 43.755 0.090 0.249 0.159 1.000

2.1.3 Use mice to synthesize with pooling rules from Drechsler

For now, we forget about the default mice imputation method, because it leads to flawed results. Therefore, we continue with the CART results only. We fit the lm model in which wgt is regressed on age and hgt.

## pool.lm.syn

## If you have a multiply imputed synthetic dataset, that is, multiple fully imputed 
## datasets that are generated by mice, using an all-ones matrix as the where matrix, 
## this function can be used to analyse the synthetic data by means of a linear 
## regression model, and to pool the results.

pool.syn <- function(mira) {
  
  if(class(mira)[1] == "mira") {
    fitlist <- mira %$% analyses
  }
  else {
    fitlist <- mira
    }
  
  m <- length(fitlist)
  
  pooled <- fitlist %>% 
    map_dfr(broom::tidy) %>%
    group_by(term) %>%
    summarise(est     = mean(estimate),
              bm      = sum((estimate - est)^2) / (m - 1),
              ubar    = mean(std.error^2),
              var_u   = (1 + 1/m) * bm - ubar,
              var     = if_else(var_u > 0, var_u, ubar),
              df      = max(1, (m - 1) * (1 - ubar / (bm + bm/m))^2),
              lower   = est - qt(.975, df) * sqrt(var),
              upper   = est + qt(.975, df) * sqrt(var), .groups = 'drop')
  pooled
}

pool2.syn <- function(mira) {
  
  if(class(mira)[1] == "mira") {
    fitlist <- mira %$% analyses
  }
  else {
    fitlist <- mira
  }
  
  m <- length(fitlist)
  
  pooled <- fitlist %>% 
    map_dfr(broom::tidy) %>%
    group_by(term) %>%
    summarise(est     = mean(estimate),
              bm      = sum((estimate - est)^2) / (m - 1),
              ubar    = mean(std.error^2),
              var_u   = (1 + 1/m) * bm - ubar,
              var     = if_else(var_u > 0, var_u, ubar),
              df      = max(m - 1, (m - 1) * (1 - ubar / (bm + bm/m))^2),
              lower   = est - qt(.975, df) * sqrt(var),
              upper   = est + qt(.975, df) * sqrt(var), .groups = 'drop')
  pooled
}

pool3.syn <- function(mira) {
  
  if(class(mira)[1] == "mira") {
    fitlist <- mira %$% analyses
  }
  else {
    fitlist <- mira
  }
  
  m <- length(fitlist)
  
  pooled <- fitlist %>% 
    map_dfr(broom::tidy) %>%
    group_by(term) %>%
    summarise(est     = mean(estimate),
              bm      = sum((estimate - est)^2) / (m - 1),
              ubar    = mean(std.error^2),
              var     = ubar + bm/m,
              df      = (m - 1) * (1 + (ubar * m)/bm),
              lower   = est - qt(.975, df) * sqrt(var),
              upper   = est + qt(.975, df) * sqrt(var), .groups = 'drop')
  pooled
}

ci_cov <- function(pooled, true_fit) {
  
  coefs <- coef(true_fit)
  vars   <- diag(vcov(true_fit))
  
  nsim <- nrow(pooled) / length(coefs)
  
  pooled %>% mutate(true_coef = rep(coefs, nsim),
                    true_var  = rep(vars, nsim),
                    cover     = lower < true_coef & true_coef < upper) %>%
    group_by(term) %>%
    summarise("True Est" = unique(true_coef),
              "Syn Est"  = mean(est),
              "Bias"     = mean(est - true_coef),
              "True SE"  = unique(sqrt(true_var)),
              "Syn SE"   = mean(sqrt(var)),
              "df"       = mean(df),
              "Lower"    = mean(lower),
              "Upper"    = mean(upper),
              "CIW"      = mean(upper - lower),
              "Coverage" = mean(cover), .groups = "drop")
}
models <- syns_cart %>% map(function(x) x %$% lm(wgt ~ age + hgt))

Then, we use the custom pooling function to pool the estimates.

pooled <- models %>% map_dfr(pool.syn)

And another custom function to summarise all results, and obtain the confidence interval coverage.

ci_cov(pooled, truemodel) %>%
  kable(digits = 3,
        caption = "Imputation results of synthesizing by means of CART, with the pooling rules by Drechsler.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Imputation results of synthesizing by means of CART, with the pooling rules by Drechsler.
term True Est Syn Est Bias True SE Syn SE df Lower Upper CIW Coverage
(Intercept) -9.216 -7.602 1.614 2.033 1.811 72.945 -20.997 5.793 26.789 1
age 2.327 2.463 0.136 0.189 0.169 91.939 1.106 3.820 2.714 1
hgt 0.191 0.169 -0.022 0.028 0.025 112.809 -0.027 0.365 0.392 1

It can be seen that the coverage of the confidence interval is nearly equal to one. Since we currently only want to make inferences with regard to the sample, this is okay. Namely, if the CI coverage with regard to the sample would be \(95\%\), the CI coverage with regard to the population would be \(0.9025\%\). Inferences with regard to the population will be discussed in the next paragraph.

However, what is questionable, is the fact that the confidence interval width has increased, even though the variance has decreased. Since the confidence interval width is a function of only the standard error and the degrees of freedom, this must be due to the fact that the degrees of freedom are smaller than the degrees of freedom as returned by the mice pooling function.

mice_df <- data.frame(term = cart_fit$var, df = cart_fit$df)
syn_df <- data.frame(term = pooled$term, df = pooled$df)

df <- bind_rows("Mice" = mice_df, "Synthetic" = syn_df, .id = "Method")

ggplot(data = df, aes(x = df, color = Method, fill = Method)) +
  geom_density(alpha = .7) +
  facet_wrap(~ term, nrow = 1, scales = "free") +
  scale_color_viridis_d(begin = .30, end = .65) +
  scale_fill_viridis_d(begin = .30, end = .65) +
  theme_classic()


2.1.4 Different synthesizing sequence

Now we have established that the mice synthesizing algorithm works as good as the synthpop synthesizing algorithm, we can check whether changing the order in which hgt and wgt are synthesized influences the result of the algorithm.

cart_wgt_hgt %>% 
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  ci_cov(truemodel) %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
term True Est Syn Est Bias True SE Syn SE df Lower Upper CIW Coverage
(Intercept) -9.216 -7.635 1.581 2.033 1.800 40.577 -20.105 4.834 24.939 0.994
age 2.327 2.459 0.132 0.189 0.165 43.528 1.171 3.747 2.576 0.996
hgt 0.191 0.170 -0.021 0.028 0.025 34.469 -0.015 0.354 0.369 0.998

2.2 Bootstrapped boys data

true_results <- bootstrap_boys %>%
  map(~ lm(wgt ~ age + hgt, .x)) %>%
  map_dfr(~ broom::tidy(.x, conf.int = TRUE)) %>%
  mutate(true_coef = rep(coef(truemodel), nsim),
         true_se   = rep(sqrt(diag(vcov(truemodel))), nsim),
         cover     = conf.low < true_coef & true_coef < conf.high) %T>% 
  assign("boot_true", ., envir = globalenv()) %>%
  group_by(term) %>%
  summarise("True Est" = unique(true_coef),
            "Boot Est"  = mean(estimate),
            "Bias"     = mean(estimate - true_coef),
            "True SE"  = unique(true_se),
            "Boot SE"   = mean(std.error),
            "DF"       = 745,
            "Lower"    = mean(conf.low),
            "Upper"    = mean(conf.high),
            "CIW"      = mean(conf.high - conf.low),
            "Coverage" = mean(cover)) %>%
  kable(digits = 3,
        caption = "Bootstrapped estimates of regression coefficients, with mean standard error and coverage based on the individual estimates of the standard error.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

synthetic_results <- boot_cart %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %T>% 
  assign("boot_syn", ., envir = globalenv()) %>%
  ci_cov(truemodel) %>%
  kable(digits = 3,
        caption = "Synthetic (CART) results with pooling rules by Drechsler.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

norm_results <- boot_cart %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  mutate(df = NA,
         lower = est - qnorm(.975) * sqrt(var),
         upper = est + qnorm(.975) * sqrt(var)) %>%
  ci_cov(truemodel) %>%
  kable(digits = 3, caption = "Normal CI approximation") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Bootstrapped estimates of regression coefficients, with mean standard error and coverage based on the individual estimates of the standard error.
term True Est Boot Est Bias True SE Boot SE DF Lower Upper CIW Coverage
(Intercept) -9.216 -9.301 -0.084 2.033 2.025 745 -13.275 -5.326 7.949 0.934
age 2.327 2.314 -0.013 0.189 0.188 745 1.945 2.683 0.738 0.932
hgt 0.191 0.193 0.001 0.028 0.028 745 0.138 0.247 0.109 0.922
Synthetic (CART) results with pooling rules by Drechsler.
term True Est Syn Est Bias True SE Syn SE df Lower Upper CIW Coverage
(Intercept) -9.216 -8.113 1.104 2.033 1.828 147.690 -19.472 3.247 22.719 0.968
age 2.327 2.418 0.091 0.189 0.167 61.161 1.289 3.547 2.258 0.964
hgt 0.191 0.176 -0.015 0.028 0.024 92.660 0.014 0.339 0.325 0.958
Normal CI approximation
term True Est Syn Est Bias True SE Syn SE df Lower Upper CIW Coverage
(Intercept) -9.216 -8.113 1.104 2.033 1.828 NA -11.696 -4.529 7.167 0.860
age 2.327 2.418 0.091 0.189 0.167 NA 2.091 2.744 0.653 0.818
hgt 0.191 0.176 -0.015 0.028 0.024 NA 0.128 0.224 0.096 0.790

We find that some bias has been introduced due to synthesizing the relationships, so we investigate this in further depth. Additionally, the coverage is somewhat too high, and the average confidence interval width has increased dramatically. Therefore, we will further investigate where this increase stems from.

plot_dat <- bind_rows("Synthetic" = boot_syn %>% 
                                    mutate(ciw = upper - lower) %>%
                                    select(term, est, ciw = ciw), 
                      "Real values" = boot_true %>% 
                                      mutate(ciw = conf.high - conf.low) %>%
                                      select(term, est = estimate, ciw = ciw),
                      .id = "Method")


ggplot(plot_dat, aes(x = est, color = Method, fill = Method)) +
  geom_density(alpha = .7) +
  facet_wrap(~ term, nrow = 1, scales = "free") +
  scale_color_viridis_d(begin = .30, end = .65) +
  scale_fill_viridis_d(begin = .30, end = .65) +
  theme_classic() +
  ggtitle("Distribution of the bootstrapped and synthetic regression coefficients.")

It can be seen that the bias is systematic, and thus not due to any extreme values that pull the average estimate to zero. Additionally, it seems that the pooled standard error indicated in the table is somewhat too low for the synthetic data. Namely, it appears in the plot that the spread in the synthetic data is roughly equal to the spread in the bootstrapped, but real, data. However, in the table above, the synthetic standard errors are much smaller than the mean over the manually calculated standard error of the bootstrapped estimates. When the standard errors are calculated over the obtained regression estimates, we find that the synthetic standard error increases much more when we calculate it over the observed synthetic estimates.

boot_se <- boot_true %>% group_by(term) %>% summarise("Mean bootstrapped SEs" = sd(estimate))
boot_syn_se <- boot_syn %>% group_by(term) %>% summarise("Mean synthetic SEs" = sd(est))

bind_cols(boot_se, boot_syn_se[,2]) %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)
term Mean bootstrapped SEs Mean synthetic SEs
(Intercept) 2.105 2.094
age 0.207 0.209
hgt 0.030 0.030
ggplot(plot_dat, aes(x = ciw, color = Method, fill = Method)) +
  geom_density(alpha = .7) +
  facet_wrap(~ term, nrow = 1, scales = "free") +
  scale_color_viridis_d(begin = .30, end = .65) +
  scale_fill_viridis_d(begin = .30, end = .65) +
  theme_classic() +
  ggtitle("Distribution of the confidence interval width (CIW) of the bootstrapped and the synthetic datasets.")

When we look at the confidence interval width, we find that the bootstrap CI is of approximately the same length in all iterations, while the CI’s of the synthetic results spread out much more.


2.3 Improving the bias/variance tradeoff

The fact that there is bias introduced, might be due to the fact that there is too little variance in the estimates of the synthetic values. Therefore, an approach that increases the variance of the estimates might yield estimates that display less bias. This can be achieved by decreasing the number of iterations, that is, setting the mice parameter maxit to 1, instead of 5. Additionally, the complexity parameter can be set to \(10^{-8}\) instead of \(10^{-4}\), so that more splits will be made. Also, the value of the rpart parameter minbucket, that is, the minimum number of observations in any terminal node, is set to 3 instead of 5, so that, once again, more splits can be made. By means of these changes, the within variance increases, while the between-imputations variability decreases. The results with 500 iterations on the same, singly imputed boys dataset can be seen below.

syns_cart_maxit_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  ci_cov(., truemodel) %>%
  kable(digits = 3,
        caption = "Synthesizing the same data.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Synthesizing the same data.
term True Est Syn Est Bias True SE Syn SE df Lower Upper CIW Coverage
(Intercept) -9.216 -8.391 0.825 2.033 1.821 52.299 -20.543 3.761 24.305 1.000
age 2.327 2.405 0.078 0.189 0.164 29.103 1.169 3.641 2.472 0.998
hgt 0.191 0.179 -0.012 0.028 0.024 38.830 0.002 0.357 0.355 1.000

Additionally, the same approach is used to synthesize the 500 bootstrapped samples from the completed boys data.

boot_cart_maxit_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  ci_cov(., truemodel) %>%
  kable(digits = 3,
        caption = "Synthesizing bootstrap samples with adjusted imputation settings.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Synthesizing bootstrap samples with adjusted imputation settings.
term True Est Syn Est Bias True SE Syn SE df Lower Upper CIW Coverage
(Intercept) -9.216 -9.080 0.137 2.033 1.952 238.879 -16.892 -1.268 15.624 0.956
age 2.327 2.332 0.006 0.189 0.178 240.207 1.484 3.180 1.696 0.952
hgt 0.191 0.190 -0.002 0.028 0.027 215.781 0.067 0.312 0.245 0.946

However, I’m not sure whether it is sensible to set the number of iterations to maxit = 1. Therefore, I additionally included an example that shows the same approach, that is, with the complexity parameter cp = 1e-08 and the parameter minbucket = 3, but with maximum number of iterations fixed at maxit = 5.

boot_cart_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  ci_cov(., truemodel) %>%
  kable(digits = 3,
        caption = "Synthesizing bootstrap samples with adjusted imputation settings but with maxit = 5.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Synthesizing bootstrap samples with adjusted imputation settings but with maxit = 5.
term True Est Syn Est Bias True SE Syn SE df Lower Upper CIW Coverage
(Intercept) -9.216 -9.037 0.180 2.033 1.964 250.209 -16.683 -1.390 15.293 0.956
age 2.327 2.337 0.010 0.189 0.178 218.191 1.513 3.160 1.647 0.940
hgt 0.191 0.189 -0.002 0.028 0.027 177.875 0.071 0.306 0.235 0.954

It can be seen that the bias indeed decreases, while the standard error of the estimates increases. This is due to the fact that the within-imputation variance increases, while the between-imputation variance decreases somewhat.

cart_default <- boot_cart %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  group_by(term) %>%
  summarise(Within = mean(ubar), 
            Between = mean(bm), 
            Total = mean(var),
            .groups = 'drop')

cart_adj_maxit <- boot_cart_maxit_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  group_by(term) %>%
  summarise(Within = mean(ubar), 
            Between = mean(bm), 
            Total = mean(var),
            .groups = 'drop')

cart_adj <- boot_cart_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool.syn) %>%
  group_by(term) %>%
  summarise(Within = mean(ubar), 
            Between = mean(bm), 
            Total = mean(var),
            .groups = 'drop')

bind_rows("Cart default" = cart_default,
          "Cart with adjusted settings and maxit = 1" = cart_adj_maxit,
          "Cart adjusted settings" = cart_adj, .id = "Method") %>%
  kable(digits = 3,
        caption = "Within, between and total imputation variance for default CART, and adjusted CART.") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Within, between and total imputation variance for default CART, and adjusted CART.
Method term Within Between Total
Cart default (Intercept) 3.890 2.164 3.516
Cart default age 0.033 0.022 0.030
Cart default hgt 0.001 0.000 0.001
Cart with adjusted settings and maxit = 1 (Intercept) 4.053 1.447 3.892
Cart with adjusted settings and maxit = 1 age 0.035 0.015 0.033
Cart with adjusted settings and maxit = 1 hgt 0.001 0.000 0.001
Cart adjusted settings (Intercept) 4.054 1.374 3.934
Cart adjusted settings age 0.035 0.014 0.033
Cart adjusted settings hgt 0.001 0.000 0.001

3 Adjusted degrees of freedom

Next to the analyses displayed previously, I adjusted the degrees of freedom to \(\text{max}(m - 1, \nu)\), as is done in Reiter and Drechsler (2010). This yields the following results for the initial analyses, that is, bootstrapping with the mice default settings, but with the method set to cart. Note that this concerns 500 bootstrap samples from the complete data.

boot_cart %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool2.syn) %T>% 
  assign("boot_syn", ., envir = globalenv()) %>%
  ci_cov(truemodel) %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
term True Est Syn Est Bias True SE Syn SE df Lower Upper CIW Coverage
(Intercept) -9.216 -8.113 1.104 2.033 1.828 149.153 -12.616 -3.610 9.006 0.906
age 2.327 2.418 0.091 0.189 0.167 62.838 1.998 2.837 0.839 0.886
hgt 0.191 0.176 -0.015 0.028 0.024 94.288 0.115 0.237 0.122 0.866

Additionally, the results can be computed for the data with the smaller bias settings. That is, setting maxit = 1, cp = 1e-08 and minbucket = 3.

boot_cart_maxit_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool2.syn) %>%
  ci_cov(., truemodel) %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
term True Est Syn Est Bias True SE Syn SE df Lower Upper CIW Coverage
(Intercept) -9.216 -9.080 0.137 2.033 1.952 239.625 -13.563 -4.597 8.966 0.940
age 2.327 2.332 0.006 0.189 0.178 241.218 1.911 2.753 0.842 0.914
hgt 0.191 0.190 -0.002 0.028 0.027 216.711 0.127 0.252 0.125 0.908

And, once again, the same analysis as the previous one, but with maxit = 5.

boot_cart_cp_min3 %>%
  map(function(x) x %$% lm(wgt ~ age + hgt)) %>%
  map_dfr(pool2.syn) %>%
  ci_cov(., truemodel) %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
term True Est Syn Est Bias True SE Syn SE df Lower Upper CIW Coverage
(Intercept) -9.216 -9.037 0.180 2.033 1.964 250.899 -13.523 -4.550 8.973 0.936
age 2.327 2.337 0.010 0.189 0.178 219.162 1.917 2.756 0.839 0.902
hgt 0.191 0.189 -0.002 0.028 0.027 178.770 0.127 0.251 0.125 0.930

References

Van Buuren, S., & Groothuis-Oudshoorn, K. (2011). mice: Multivariate imputation by chained equations in r. Journal of Statistical Software, 45(3), 1–67. Retrieved from https://www.jstatsoft.org/v45/i03/